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When using black hole excision to numerically evolve a black hole spacetime with 
no continuous symmetries, most 3 + 1 finite differencing codes use a Cartesian grid. 
It's difficult to do excision on such a grid because the natural r — constant excision 
surface must be approximated either by a very different shape such as a contained 
cube, or by an irregular and non-smooth "LEGO^ sphere" which may introduce 
numerical instabilities into the evolution. In this paper I describe an alternate 
scheme, which uses multiple {r x (angular coordinates)} grid patches, each patch 
using a different (nonsingular) choice of angular coordinates. This allows excision 
on a smooth r = constant 2-sphere. 

I discuss the key design choices in such a multiple-patch scheme, including the choice 
of ghost-zone versus internal-boundary treatment of the interpatch boundaries (I 
use a ghost-zone scheme), the number and shape of the patches (I use a 6-patch 
"inflated-cube" scheme), the details of how the ghost zones are "synchronized" by 
interpolation from neighboring patches, the tensor basis for the Einstein equations 
in each patch, and the handling of non-tensor field variables such as the BSSN 
(I use a scheme which requires ghost zones which are twice as wide for the BSSN 
conformal factor as for and the other BSSN field variables). 
I present sample numerical results from a prototype implementation of this scheme. 
This code simulates the time evolution of the (asymptotically fiat) spacetime around 
a single (excised) black hole, using 4th order finite differencing in space and time. 
Using Kerr initial data with J/rri^ = 0.6, I present evolutions to t ^ 1500m. 
The lifetime of these evolutions appears to be limited only by outer boundary 
instabilities, not by any excision instabilities or by any problems inherent to the 
multiple-patch scheme. 
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I. INTRODUCTION 

When time-evolving a black hole spacetime in 3 -f 1 (Cauchy) numerical relativity, the 
numerical computation must somehow be kept from encountering the singularity(ies) within 
the black hole(s). There are three common means of doing this: freezing slicings (Lichnerow- 
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icz (1944); Smarr and York (1978)), the Brandt-Briigmann "puncture" method (Brandt and 
Briigmann (1997)), and excision (Anninos et al. (1995); Seidel and Suen (1992); Thorn- 
burg (1985, 1987, 1993); Unruh (1984)). In this paper I focus on excision in the context 
of finite-difference evolutions, in particular the case where the spacetime has no continuous 
symmetries and hence a fully 3-spatial-dimensional numerical grid must be used. For con- 
venience of exposition I primarily consider the case where there is precisely one (excised) 
black hole in each slice; I briefly discuss possible extensions to multiple black holes in the 
conclusions. 

Polar spherical coordinates and grids centered on a black hole provide a natural way to 
do excision, with an r = constant excision surface being a (smooth) grid plane. However, 
problems with z axis coordinate singularities (and the difficulty of generalizing to multiple 
black holes) have led most researchers to switch to Cartesian grids over the past decade. 

Unfortunately, it's difficult to do excision on a Cartesian grid: 

1. If the excision boundary is a cube, there are severe causality problems (Calabrese and 
Neilsen (2004b); Lehner (2003); Lehner et al. (2004a); Scheel and Kidder (2000)). For 
Schwarzschild spacetime in Eddington-Finkelstein (spacehke) coordinates, the excised 
region must be very small to keep the evolution causal. For Kerr spacetime in Kcrr- 
Schild coordinates this problem is even worse, and for most spin parameters there is 
no cubical excision region which keeps the evolution causal. For a dynamic spacetime 
this problem is even more severe. 

2. If the excision boundary is some other regular polyhedron which better approximates 
a sphere (such as the 14-sided shape obtained by cutting a pyramid off each of a cube's 
corners) , the continuum causality problem is less severe, but the additional edges and 
corners considerably complicate the construction of a stable finite differencing scheme. 

3. If the excision boundary is an irregular "LEGO sphere" , then at the continuum level 
(approximating the boundary as a sphere centered on the origin), any non-extremal 
Kerr spacetime admits a causal evolution (Thornburg (1993, figure 3.8)).^'^ However, 
at the discrete level, the irregularity of the boundary makes it very difficult to con- 
struct stable finite differencing schemes. (Notably, it's not clear whether the powerful 
numerical schemes of Calabrese et al. (2003a,b) can be extended to this case.) The 
causality problem of case 1 may also still apply to the individual "LEGO sphere" grid 
points on the excision boundary. 

In this paper I discuss an alternative technique for doing excision, which is similar to the 
use of a polar spherical grid, but with the {9^ (p) angular grid replaced by multiple angular 
grid patches which collectively cover S'^ in a non-singular manner. The corresponding {r x 
(multiple angular patches)} coordinate system can be constructed to have no coordinate 
singularities, and the multiple-patch grid still has r = constant as a (smooth) excision 
surface. 



^ For the Kerr slicing of Kerr spacetime, an r = constant excision boundary gives a causal evolution 
whenever the boundary lies between the inner and outer horizons. This gives a relatively broad range of 
allowable boundary positions for small spin parameters, narrowing as the spin increases. 

^ Since the causality depends only algebraically on the 4-metric, the same property also holds for all "nearby" 
spacetimes in some functional neighborhood of Kerr spacetime. 
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Like a polar spherical grid, such a multiple-patch grid provides a smooth r = constant 
outer boundary for the imposition of outgoing radiation conditions, and makes it easy to 
place the r = constant shells of grid points nonuniformly in radius so as to give high 
resolution close to the black hole while still having the outer boundary relatively far away.^ 

Multiple-patch finite difference schemes with smooth excision surfaces have long been 
used in computational fluid dynamics (see, for example, the discussions of Rubbert and 
Lee (1982), Thompson et al. (1985, section IL4), Chesshire and Henshaw (1990, 1994), 
Gustafsson et al. (1995, section 13.4), Brown et al. (1997), and Petersson (1999)). However, 
such schemes have seen relatively little use in numerical relativity, and that mainly for 
elliptic problems, either initial data construction (Thornburg (1985, 1987)) or apparent- 
horizon finding (Thornburg (2004)).*^ 

Bishop et al. (1996) have described a multiple-patch scheme using a pair of overlap- 
ping stereographic-coordinate patches to cover S"^, with field variables in each patch's ghost 
zone interpolated (in 2-D) from the other patch. Their original application was Cauchy- 
characteristic matching for radiation boundary conditions, but the same technique has since 
been applied for fully-nonlinear null-cone evolutions (see, for example. Bishop et al. (1997); 
Gomez et al. (1998)). 

Calabrese and Neilsen (2004a,b) have described a multiple-patch scheme for the axisym- 

metric evolution of a scalar field on a Schwarzschild background. They use two overlapping 
patches: a polar spherical inner patch (with a smooth r = constant excision surface) and a 
cylindrical outer patch. Ghost-zone values of all the dynamical fields in each patch's ghost 
zone are set by (2-D) bilinear interpolation from the other patch. Their discretization (Cal- 
abrese et al. (2003a,b)) is based on finite difference operators which preserve discrete energy 
norms within each patch, maximally dissipative boundary conditions, and added artificial 
dissipation to ensure stability. Their scheme gives long-term stable evolutions with 2nd order 
overall (global) accuracy;^ a variant discretization gives 3rd order overall accuracy. 

Lehner et al. (2004b); Reula (2003); Tigho (2003, 2004) have described a different 
multiple-patch scheme for studying scalar field tails in Schwarzschild or Kerr spacetime. 
They use an "internal boundaries" scheme, where adjacent patches just touch (rather than 
overlapping). By using a symmetric hyperbolic form of the PDEs, the dynamical fields can 
be unambiguously decomposed into ingoing and outgoing modes crossing each interpatch 
boundary. By adding penalty terms to the field equations, the difference between the value 
of each ingoing mode at a patch boundary and its value as an outgoing mode in the adjacent 
patch, can be driven to zero. 



^ Berger-Oliger mesh refinement (Berger (1982, 1986); Berger and Colella (1989); Berger and Oliger (1984)) 
provides a much more general and efficient solution here, but it's quite difficult to implement for the 
Einstein equations (see, for example, Briigmann et al. (2004); Choptuik (1986, 1989); Imbiriba et al. 
(2004); Licbling (2002); Schnetter et al. (2004), and references therein). 

* Multiple-patch spectral methods have been used in numerical relativity by a number of researchers (see, 
for example, Ansorg et al. (2004, 2003); Grandclement et al. (2001, 2002); Kidder and Finn (2000); Pfeiffer 
et al. (2003), and references therein). 

^ Although bilinear interpolation has an 0((Ax)^) tnmcation error, this error is non-smooth at the inter- 
polation points (Thornburg (1999a, appendix F)), so spatial derivatives of the field variables are generi- 
cally only 1st order accurate at grid points where interpolated values enter into the derivative molecules. 
However, in practice this doesn't seem to seriously affect the overall (global) convergence of the scheme 
(Calabrese and Neilsen (2004a)). I discuss this point in more detail in section III. 
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In this paper I describe a multiple-patch excision scheme for the full nonlinear vacuum 
3 + 1 Einstein equations. In terms of past multiple-patch schemes in numerical relativity, 
this scheme most resembles those of Gomez et al. (1997) and Calabrese and Neilsen (2004b), 
in that I couple the patches by ghost-zone interpolation of all the dynamical fields. However, 
my scheme uses patches with may either overlap or just touch, and there are many other 
differences in detail. 

I first implemented the ^^g-K" ADM form of the 3-1-1 Einstein equations (York (1979)), 
but it's now well known that these are ill-posed even in the absence of boundary conditions 
(see, for example, the numerical tests of Alcubierre et al. (2004)), and I saw rapid error 
growth near the interpatch boundaries, with evolutions crashing in ^ 100m (Thornburg 
(2003)). I now use the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) form of the 3 + 1 
Einstein equations first introduced by Nakamura and Oohara (1989, 1998); Nakamura et al. 
(1987); Shibata (1999); Shibata and Nakamura (1995), and popularized by Baumgarte and 
Shapiro (1999). In particular, I use the "actively forced A" variant described by Alcubierre 
et al. (2000). 

In this paper I also present initial numerical results for the BSSN system, evolving octant- 
symmetry Kerr spacetime with high accuracy for ^ 1000m, and lasting for t ^ 1500m before 
crashing. These evolutions are unstable at the outer boundary (where I have only partially 
implemented Sommerfeld outgoing-radiation boundary conditions), but they show no signs 
of any instabilities caused by the multiple-patch scheme. There are also no signs of any 
instabilities at the inner (excision) boundary. 

A. Notation 

I generally follow the sign and notation conventions of Wald (1984). I use the Penrose 
abstract-index notation, with indices a-z running over both Cartesian coordinates = 
{x, y, z) in a spacelike 3 + 1 slice, and the generic patch coordinates a;* = (r, p, a) defined in 
section II. B. 

Qij is the 3-metric in the slice, g is its determinant, and Vj is the associated covariant 
derivative operator. Kij is the extrinsic curvature of the slice (I use the sign convention of 
York (1979), not that of Wald (1984)) and K = is its trace. "7V-D" abbreviates "TV- 
dimensional" , and specifically refers only to spatial dimensions (i.e. the temporal dimension 
is never included). 

I use sans-serif font for the names of patches, ghost zones, and (in appendix A) ghost 
zone widths and patch interpolation regions. In section II. E I use indices abc and ijk for 
the (r, p, a) coordinates of a pair of neighboring patches p and q respectively, and I use the 
notation /(p) for the grid function / in the local coordinate basis of patch p. 

In the context of a particular angular patch boundary, I use angular coordinates (_L, ||) 
which are (perpendicular, parallel) to the boundary. In appendix A I also use the corre- 
sponding integer grid-point coordinates and grid-function array indices (iperp, ipar), and 
I use a notation inspired by the C++ programming language, where (p. iperp, p. ipar) refers 
to the (iperp, ipar) coordinates of patch p. 

In appendix C I use overbars (as in 0) to refer to field variables in (strictly speaking, 
their coordinate components with respect to) the Cartesian coordinates. 
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II. DESIGN OF THE MULTIPLE-PATCH SYSTEM 

In this section I discuss a variety of design issues which arise when trying to design and 
implement such an {r x (multiple angular patches)} numerical scheme for the 3 + 1 vacuum 
Einstein equations. 

A. Treatment of Interpatch Boundaries 

A fundamental design question for any multiple-patch scheme is how the interpatch 
boundaries should be treated in the numerical scheme, or equivalently, how the patches 
should be coupled together. 

In terms of finite differencing, the usual ghost-zone scheme is conceptually simple: each 
patch's nominal grid is surrounded by a molecule-radius "ghost zone" of extra grid points 
on each side (face) of the patch. Values of all the dynamical fields in the ghost zones are 
computed ("synchronized") by interpolation from the adjacent patches. Finite differencing 
can then be done at all the nominal grid points without further concern for the interpatch 
boundaries. This scheme can be used with essentially any formulation of the Einstein equa- 
tions, and also generahzes to elliptic PDEs in a reasonably straightforward and efficient 
manner (Thornburg (1985, 1987, 2004)). The main problem with the ghost-zone scheme is 
that there's relatively little analytical theory regarding its stability, although Olsson and 
Petersson (1996) have obtained some promising results for 1-D model problems, and Star- 
ius (1980) has proved stability for a 1-D model problem (as well as performed numerical 
experiments using the 2-D shallow- water equations). 

In contrast, the "internal-boundaries" scheme used by Lehner et al. (2004b); Reula (2003); 
Tiglio (2003, 2004) requires a formulation of the Einstein equations in which the modes cross- 
ing each internal boundary can be cleanly separated. This must be done everywhere on each 
interpatch boundary, including (in general) in strong-field regions. Although this scheme 
does generalize to elliptic PDEs, this generalization is quite inefficient in multiple spatial 
dimensions.^ The main advantage of the internal-boundaries scheme is its amenability to 
mathematical analysis at both the continuum and the finite differencing levels (sec, for ex- 
ample. Carpenter et al. (1999)). Also, potentially-destabilizing feedback loops from one 
patch to another and back again only occur due to nonlinear effects coupling the oppositely- 
directed modes, whereas in the ghost-zone scheme all the field variables participate directly 
in such feedback loops. 

Because of its simplicity, and for historical reasons, I have chosen the ghost-zone scheme 
for this work. 



^ Briefly, in this scheme nonhnear elhptic PDEs are first Unearized in the Newton-Kantorovich sense (Boyd 
(2000, appendix C)). Linear elliptic PDEs are then solved by introducing homogeneous solutions where all 
the interpatch field values vanish, and one particular solution for each interpatch boundary point. Thus in 
a single spatial dimension, the number of particular solutions is just twice the number of patches. But in 
multiple spatial dimensions, it's twice the number of interpatch-boundary grid points, which is generally 
quite large. 
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B. Choice of Angular Patches and Coordinates 

There are a number of ways to cover S'^ with nonsingular coordinate patches, of which 
two seem particularly interesting: 

One possibility (Bishop et al. (1996)) is to use a pair of stcrcographic-coordinatc patches, 
one covering the northern hemisphere and one the southern. This has the advantage of 
relative simplicity (there is only a single interpatch boundary), and the eth formalism (Gomez 
et al. (1997)) provides an elegant way to represent and manipulate tensor fields on S"^. 

However, these coordinates have relatively large coordinate distortion near the equator, 
which makes finite differencing somewhat less accurate there. Another potential problem 
is the coordinate choice on each patch: If Cartesian stereographic coordinates are used on 
each patch, then the grids overlap irregularly near the equator (requiring 2-D interpatch 
interpolation there), while if polar stereographic coordinates are used on each patch, then 
there are coordinate singularities at the north and south poles. 

Another possibility is to use 6 patches, covering neighborhoods of the ±x, iy, and ±2; axes 
respectively. This has the disadvantage of relative complexity: there are many interpatch 
boundaries, and there are also corners where 3 patches meet. However, this system has 
the advantage of relatively low coordinate distortion, yielding accurate finite differencing. 
Also, if the coordinates are suitably chosen (1 describe this in detail below), it's possible to 
have adjacent patches always share a common angular coordinate, so only 1-D interpatch 
interpolation is needed. 

Because of the lower coordinate distortion, and the simplicity of only needing 1-D in- 
terpatch interpolation, I use a 6-patch system here. In more detail, 1 use the "inffated- 
cube" angular 6-patch system described in Thornburg (2004): Given Cartesian coordinates 
(x, z), with the excised (black hole) region at the origin, I define 3 angular coordinates on 
S'^ based on rotation angles about the xyz coordinate axes: 

/i = rotation angle about the x axis = arctan(y/z) 

u = rotation angle about the y axis = arctan(2;/^) (1) 

(p = rotation angle about the z axis = aictan{y / x) 

where all the arctangents are 4-quadrant based on the signs of x, y, and z. 1 then define 
6 coordinate patches covering neighborhoods of the ±z, ±x, and ±y axes, using the generic 
(angular) patch coordinates (p, a) defined by 

±z patch has (p, a) = (p, u) 

±x patch has (p, a) = {u, ip) (2) 
±y patch has (p, a) — (p, (p) 

Notice that each patch's (p, a) coordinates are nonsingular throughout a neighborhood 
of the patch, and that adjacent patches always share the (common) angular coordinate 
perpendicular to their mutual boundary. The name "infiated-cube" comes from another 
way to visualize these patches and coordinates: Imagine an xyz cube with xyz grid lines 
painted on its face. Now imagine the cube to be flexible, and inflate it like a balloon, so it 
becomes spherical in shape. The resulting coordinate lines will closely resemble those for 
(p, z/, (p) coordinates. 

This set of 6 patches covers S"^ without coordinate singularities. Alternatively, if the 
spacetime has z —z reflection symmetry about the coordinate origin, then the 5 patches 



7 




FIG. 1 This figure shows a multiple-grid-patch system covering the (-I-, -|-, -|-) octant of S*^ with 
3 patches, at an angular resolution of Apcj = 5°. The +z, +x, and +y patches are shown in red, 
green, and blue respectively. The patch's nominal grids (shown in thick lines) just touch; their 
ghost zones (which are 2 points wide, and are shown in thin lines) overlap. 

+z, ±x, and ±y cover the +z hemisphere of S^. Similarly, suitable sets of 4 or 3 patches 
may be used to cover quadrants or octants. Figure 1 shows an example of a 3-patch system 
covering the (+, +, +) octant of 5^. 

Defining the usual radial coordinate r = (x^ + + z'^Y^'^, the corresponding 6-patch (or 
5-, 4-, or 3-patch) {r,p,a) coordinates cover all of 3?'^ (or the corresponding subsets) with 
no singularities except at the origin (which is excised). 

C. Choice of Tensor Basis 

For numerical solution I decompose the 3-1-1 Einstein equations into coordinate compo- 
nents with respect to a tensor basis. This raises the question of what basis should be used 
in each patch. There are two natural approaches: 

Use the same (typically Cartesian) basis in each patch 

This approach avoids the need for any change of basis in the interpatch interpolation. 
By using a Cartesian basis, it may also make it easier to interface with other numerical 
relativity software which uses this type of tensor basis. 

However, this approach has the disadvantage that a nontrivial (frame) transformation 
is required at each grid point to convert grid finite differences into approximate coor- 
dinate partial derivatives. This may both slow down the code, and possibly introduce 
time evolution instabilities, or at least complicate stability analysis. It also makes it 
difficult to reuse existing software modules for the evolution. 
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Using a Cartesian basis also makes it harder to treat the radial coordinate specially for 
such things as special finite differencing near the excision boundary, outer boundary 
conditions, etc. 

Use each patch's own (r, p, cr) coordinates to define its tensor basis 

Because this approach uses a different tensor basis in each patch, it requires a change 
of basis in the interpatch interpolation. (This is discussed in detail in section II. E.) 
However, this extra computation is only needed at the patch boundaries - in the 
interior of each patch there is no extra overhead, and existing unigrid software can 
potentially be reused with little or no change. 

This approach also makes it easy to treat the radial coordinate specially in finite 
differencing and/or in the continuum formulation of the equations. 

In this work I use the second approach: in each patch I use that patch's own (r, p, a) 
coordinates to define the tensor basis. 

D. Synchronizing Ghost Zones 

As described in section II. A, I use the usual "ghost zone" technique for handhng finite 

differencing near the patch boundaries. I refer to the non-ghost-zone part of a patch's grid 
as its "nominal" grid, and to the process of computing values for a set of grid functions in 
all ghost zones of all patches as "synchronizing" these grid functions. This must be done for 
all grid functions to which finite difference molecules will be applied, i.e. in practice, for all 
grid functions representing dynamical field variables. 

Each patch is a rectangular solid ("cuboid") in its own (r, p, a) coordinates, so it has 
6 ghost zones: 2 radial (inner and outer) and 4 angular (pmin, Pmax, c"mm, and cTmax)- There 
are three types of ghost zones, each with corresponding synchronization techniques: 

A radial ghost zone is synchronized by extrapolation from the nominal grid. I describe 
this in detail in appendix A. 

A "symmetry" (angular) ghost zone is one where the patch boundary is a discrete 
symmetry plane of the spacetime, such as the x = 0, y = 0, or z = plane in the 
octant-symmetry example of figure 1. In this case the symmetry operation defines a 
mapping from the ghost zone into the nominal grid of some (possibly other) patch, so 
the ghost-zone grid function values can be computed by copying from the symmetry- 
image grid points. 

An "interpatch" (angulcir) ghost zone is one where the ghost zone lies within the nom- 
inal grid of some other patch. ^ In this case the ghost-zone grid function values can be 
computed by interpolating from the neighboring patch. Because all patches share a 
common radial coordinate, and (by construction) adjacent patches share the angular 
coordinate perpendicular to their mutual boundary, the (1-D) interpatch interpolation 
is done (independently in each line of constant perpendicular coordinate and r) only 
in the direction parallel to the boundary. 



To ensure that this is always so, adjacent patches' nominal grids must either overlap, just touch, or be 
separated by a gap of at most one grid spacing. 
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Because off-centcrcd interpolations are relatively inaccurate, and have large phase errors 
which may lead to finite differencing instabilities in the time evolution, I try to keep the 
interpatch interpolations as centered as possible. 

A major complication in ghost zone synchronization is that ghost zones have corners, i.e. 
there are ghost-zone points which are outside the nominal grid in more than one dimension. 
It's quite tricky to ensure that the corner ghost-zone points are properly updated while still 
keeping the interpatch interpolations as centered as possible near the patch corners. I use the 
3- phase algorithm of Thornburg (2004, appendix A), generahzed to also handle radial ghost 
zones and non-scalar grid functions, to synchronize ghost zones. I describe the resulting 
algorithm in detail in appendix A. It turns out that if the interpatch-interpolation molecule 
is 4-point or smaller, then all the interpatch interpolations can in fact be kept centered. If the 
interpatch-interpolation molecule is 5-point or larger, then all the interpatch interpolations 
can still be kept centered except in the immediate neighborhood of a radial line where 
3 patches meet. 



E. Change of Tensor Basis 

In numerical relativity we need to deal with tensors and other non-scalar grid functions, 
so interpatch symmetry operations and interpolations entail a change of basis (step 4 in the 
synchronization algorithm of figure 9). 

Without loss of generality I assume that all the fields are known at a point (event) in 
some patch p in some other patch q's (r, p, cr) basis {a;*(q)}, and we wish to transform them 
to patch p's (r, p, cr) basis {a;"(p)}. I define the transformation matrices 

" ai(p)"ai(p)' ^ ' 

These are straightforward to compute analytically from the coordinate definitions (1). 
Of the BSSN field variables, a and K are scalars, /3* is a tensor and transforms as 

/3'^(p)=X\/3^(q) (4) 

and (p is (proportional to) the logarithm of a tensor density and transforms as 

0(p)=0(q) + ilog|r"||| (5) 

where || refers to each patch's angular coordinate parallel to p and q's mutual interpatch 
boundary, gij and Aij are tensor densities, and transform as 

~gM = |r"|ir'^V^r^^,,-(q) (6a) 
AM = \Y\\\~'''''Y\Y\Ai^{q) (6b) 



10 



f * isn't a tensor or tensor density, so it transforms in a more complicated manner (I 
outline the derivation of this in appendix B), 

^(p) = \Y\^\'^' X\f'^{q) + X\Y\J''ip) 

- 2 |Fll||f/'xV'(q)5^0(q) + 2r\p)d,cP{p) (7) 

The first two terms here are straightforward to implement, but the 3rd and 4th terms are 
problematic in several ways: 

• The 3rd term of (7) involves the inverse conformal metric ^'^^(q). Since this is a 
patch-q quantity, it's presumably only known on the patch-q grid, and so must be 
interpolated to the (incomeasurate) patch- p grid. This requires either keeping g^^ as a 
grid function so that it can be interpolated (instead of the more economical-of-memory 
alternative of only storing its value at the current grid point), or interpolating Qij and 
then computing g^^ from that. In my code I do the latter. 

• The 3rd term of (7) also involves 5^0(q). Since is a dynamic field variable, this partial 
derivative must be computed by finite differencing. This is difficult because this is a 
patch-q derivative, but it's needed at a patch-p grid point. The computation can be 
done by either using a differentiating interpolator*, or (less efficiently) by storing the 
di4>{c{) as grid functions and interpolating them directly. (This computation is being 
done at a point which lies within patch q's nominal grid, so there's no problem in 
computing ^^^(q).) For historical reasons, I use the latter technique in my current 
code. 

• The 4th term of (7) involves ^^^(p), which is particularly difficult to compute, because 
(as in the 3rd term) this derivative must be computed by finite differencing, and this 
must be done at a point which lies in patch p's ghost zone. There are two ways to 
compute this term: 

— i9f,0(p) can be computed by straightforward finite differencing of 0(p) if has a 
ghost zone twice as wide as that of and the other BSSN field variables, so that 
the computation point (in F^'s ghost zone) is still guaranteed to be surrounded 
by at least a molecule-sized neighborhood of 0(p) grid points. 

— Alternatively, using the analytically- known 0(p) <-> 0(q) transformation (5), 
db(f){p) can be rewritten in terms of the d(d>{c\), which can be computed using 
the methods described above for the 3rd term of (7).^ This avoids having to 
introduce double-width ghost zones for (p. 



An interpolator generally works by (conceptually) locally fitting a fitting function (usually a low-degree 
polynomial) to the data points in a neighborhood of the interpolation point, then evaluating the fitting 
function at the interpolation point. By evaluating the derivative of the fitting function, the 9^(/'(q) values 
can be interpolated very cheaply, using only the input (/)(q) values which are used anyway for interpolating 
^(q) to the patch-p grid. 

Actually this is only needed for the perpendicular derivative d±(f){p), since the radial and parallel deriva- 
tives dr(t){p) and 9||(p)0(p) are taken in directions parallel to the ghost zone, and so can be computed by 
standard finite differencing techniques without using any ^(p) data from outside the ghost zone. 
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Because it minimizes the amount of interpolation to be done, and seems to give a 
slightly simpler and more easily understood numerical scheme, I use the double-width 
ghost zone approach for all the numerical results presented here. In particular (given 
the 5-point finite difference molecules described in section III), I use angular ghost 
zones which are 2 grid points wide for most grid functions, and 4 grid points wide 
for (f). 

Although this approach is conceptually simple, it turns out to be somewhat awkward 
to retrofit to an existing code, as large parts of the code generally implicitly assume 
that the ghost-zone width is an inherent property of a grid, rather than (potentially) 
varying from one grid function to another. 

III. FINITE DIFFERENCING 

In this section I describe the finite differencing scheme used to obtain the numerical 
results presented in section IV. This is a straightforward generalization to 3-D of the 1-D 
finite differencing scheme I described in Thornburg (1999b). 

To allow higher resolution near the black hole while allowing the outer boundary to be 
placed relatively far away, I place the r = constant shells of grid points non-uniformly in 
r: they are uniformly spaced in a new (dimensionless) radial coordinate w = w{r). The 
implementation of this nonuniform gridding is identical to that of Thornburg (1999b), with 
the parameters tq = 1.5m, a = oo, b — 5m, and c = 100m. For these parameters w 
qualitatively resembles a logarithmic radial coordinate in the inner part of the grid, and a 
uniform radial coordinate in the outer part of the grid. 

I use the method of lines, with a low-storage variant of the classical 4th order Runge-Kutta 
time integrator (Blum (1962); Williamson (1980)). For my nonuniform-grid parameters and 
coordinate conditions (section IV. B), the empirical CFL limit of my code is a Courant 
number of At/ Aw = 0.63 ± 0.01; I use At/ Aw — 0.5 for all the numerical results presented 
here. 

For the spatial finite differencing, I first transform all (r, p, a) coordinate partial deriva- 
tives into (w,p, cr) coordinate partial derivatives (for example, drf = {dw /dr)dwf) for any 
field variable /). I then approximate the {w,p,(j) derivatives by the usual centered 5-point 
4th order 1st or 2nd derivative molecules as appropriate, except that for the shift vector 
advection derivatives I use off-centered 5-point molecules, upwinded by 1 grid point in the 
radial direction based on the sign of (which is always positive for the numerical results 
presented here). 

At the inner and outer grid boundaries I use 4th order Lagrange polynomial extrapolation 
of the nominal-grid field variables to fill in values in the radial ghost zones (steps 1 and 6 in 
the ghost-zone synchronization algorithm of figure 9). Because of the upwind derivatives in 
the radial direction, the ghost zones are 2(3) grid points wide at the inner(outer) boundary 
for most grid functions, or 4(6) grid points wide for 0. 

I use 5th order (6-point) Lagrange polynomial interpolation for all the interpatch inter- 
polations (step 3 in the ghost-zone synchronization algorithm of figure 9). This gives an 
0((Ax)^) truncation error for all the interpolated field variables. However, these errors - 
and thus the interpolated field variables themselves - are non-smooth at the interpolation 
points (Thornburg (1999a, appendix F)). Because of this non-smoothness, the consistency 
argument of Choptuik (1991) does not apply here: taking numerical derivatives of the in- 
terpolated field variables does locally lower the order of accuracy. 
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For the ADM equations, taking numerical 2nd derivatives of the interpolated field vari- 
ables gives 4th order accuracy, and my numerical tests (omitted here in the interests of 
brevity) confirm that the scheme attains this. However, for the BSSN equations the F* in- 
terpatch change-of-basis transformation (7) involves numerical 1st derivatives of the BSSN 
conformal factor 0, so the results of this transformation have non-smooth 5th order errors 
(as well as the usual smooth 4th order errors). Taking numerical 2nd derivatives of F* in 
the BSSN evolution and constraint equations thus lowers the local accuracy of my scheme 
to 3rd order near interpatch boundaries. ^° 

The numerical results presented in section IV. E show that despite the lower local order of 
accuracy near the interpatch boundaries, in practice my scheme remains globally 4th order 
accurate everywhere away from the grid and patch boundaries. Like the analogous result 
of Calabrese and Neilsen (2004a,b) (described in footnote 5), and those of Imbiriba et al. 
(2004), this is in accordance with theoretical arguments (Gustafsson et al. (1995, page 571)) 
that under suitable conditions boundary conditions can be approximated one order lower 
in accuracy than interior equations, without affecting the global order of accuracy of the 
scheme. 

Because the inner (excision) grid boundary (placed at r = 1.5m for all the numerical re- 
sults presented here) is spacelike, no further boundary conditions are needed there, and I use 
the usual interior evolution equations to determine the time evolution of the field variables. 
Because 1 only use 4th order radial extrapolation, not 5th order, the finite differencing is 
only 3rd order accurate for 2nd derivatives at the inner and outer grid boundaries. This 
doesn't degrade the overall 4th order evolution accuracy at the inner boundary because the 
grid points there are an "outfiow" boundary with respect to the evolution causality (Gary 
(1975); Gustafsson (1971, 1975, 1982); Gustafsson and Kreiss (1979)). 

At the outer boundary I use the outer boundary conditions described in appendix G 
to determine the time derivatives of the BSSN field variables. Unfortunately, there's no 
reason to think these boundary conditions are either well-posed or constraint-preserving, 
and outer boundary instabilities seem to be the limiting factor for the evolutions I present 
in section IV. 

Olsson and Petersson (1996) found that some artificial dissipation was necessary to obtain 
stable evolutions with a two-patch scheme for a simple (linear) 1-D model problem. I 
haven't included any artificial dissipation in my numerical scheme, though a small amount of 
dissipation is present inherently due to the time integrator and the upwind finite differencing 
of the shift vector advection terms. I don't know whether adding artificial dissipation would 
improve my scheme's stability. 



This could be raised to 4th order by using a higher-order or smoother (Hermite or spUne) interpatch 
interpolation, but I have only made limited trials of this thus far. 

If it were desirable, there would be no particular problem with using 5th order radial extrapolation to 
avoid the lower order accuracy at the radial boundaries (Thornburg (1999b)); I have not done this in the 
present work only for historical reasons. 

The errors at the outer boundary are dominated by those from the continuum outer boundary conditions 
(appendix C), so the 3rd order finite differencing there is unimportant. 
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IV. NUMERICAL RESULTS 

In this section I present various sample results from a prototype implementation of the 
numerical scheme described in this paper. This code is a standalone uniprocessor code, not 
using any larger relativity toolkit such as Cactus (Goodale et al. (2003)). The multiple- 
patch infrastructure (roughly lOK lines of C++) was quite difficult to design and debug, 
though I suspect a reimplementation would be somewhat simpler. I found no problems 
in combining the multiple-patch scheme with "relativity" code machine-generated from a 
higher-level tensor form. 

A. Initial Data 

For the numerical results presented here, I use initial data for the lapse, shift, and BSSN 
field variables which is a Kerr-coordinate slice of Kerr spacetime, with spin J/m^ = 0.6. I 
use this same slice as the background slice for my outer boundary conditions (appendix C). 
In Kerr coordinates the horizon is at r = (1 + yl^^a?)m = 1.8m for this spin; as mentioned 
above, I place the inner (excision) grid boundary at r = 1.5m. 

B. Coordinate Conditions 

I use a time-independent shift vector for all the numerical results presented here, with /9* 
set to a suitable (position-dependent) value on the initial slice and not updated thereafter. 

I use generalizations of the Bona-Masso slicings (Bona et al. (1995)). In particular, I 
use a slicing condition slightly adapted from the "X-driver" condition of Alcubierre et al. 
(2003a,b), 

dta = -af{a) {aK - (8) 

with f{a) — Aa^. As discussed by Alcubierre et al. (2003a, section III. A), these slicings 
have a gauge propagation speed of 

ciapse = a^J f{a) (9) 

1 made some evolutions with A = 2, n = —1 (the "1 -|- log" slicing recommended by 
Alcubierre et al. (2003a)), but these suffered from severe gauge instabilities. 1 use A = 2, 
n = (a variant of harmonic slicing) for all the results reported here. 

C. Numerical Grid Parameters 

I have evolved this initial data with a number of different finite differencing grids, as 
shown in table I. All the results presented here are for data with octant symmetry. 

D. Diagnostics 

I use two main diagnostics to assess the code's accuracy: The first is the energy constraint 
C = R — KijK'^^ + . To determine a scale for nonzero values of this, it's useful to also 
consider Cabs = 1-^1 + l-^ij-^*"' ! + its value Cabs.exact computed analytically for my 

(Kerr) initial data. Cabs,exact is 0(1) near the black hole, but falls off ~ at large r. The 
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Radial Resolution Ar 


Position 


Resolution 


Patch 








Model 


Aw 


Inner 


w = 0.12 


Outer 




^'max 


Apa 


Overlap Outcome 






33k-wrmax4 
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TABLE I This table shows the grid parameters and run outcomes for the various evolutions. The 
radial resolution is given first as the radial grid spacing in w, then as the radial grid resolution 
Ar at the inner boundary, the w = 0.12 (r = 2.19m) position used for the angular convergence 
plots, and the outer boundary. The "patch overlap" is specified in one of two ways: ±m (n) or 
±x° (n). ibm or ±x° gives the distance of the patches' nominal grid boundaries from 45° (either 
m angular grid spacings, or x degrees), while n gives the number of perpendicular coordinate values 
common to a pair of overlapping patches. For example, for a model with angular grid resolution 
Apa = 3°, "±3 (7)" and "9° (7)" both mean that the patches' nominal grid boundaries arc at 
45° zb 3 Apa = 45° ±9°, so adjacent patches have 7 grid points in common in their perpendicular- 
to-the-boundary direction. (Thus the 50k-overlap3 and 50k-9overlap models are actually identical.) 

"relative energy constraint" C/Cabs,exact is a dimensionless measure of the accuracy with 
which the code is approximating a solution of the Einstein equations. 

My second diagnostic is the error of the BSSN state vector with respect to its analytical 
(Kerr) value for my initial data, 

SS ^ {S<l>r + Y.{5~g,,f + {5Kf + Y.{5A,f + Y,{5rf (10) 

i<j i<j i 

where for each BSSN field variable / e {0, Qij, K, A^j, F*}, Sf = f — /xerr- 

I consider a run to "crash" if any BSSN field variable exceeds 10^° in magnitude at any 
grid point. Table I lists the times at which this occurs for each run. 

E. Convergence 

As discussed in section III, in the limit of infinite resolution and in the absence of outer- 
boundary effects, my numerical scheme should give 4th order convergence for the energy 
constraint C in patch interiors, 3rd order for C near interpatch boundaries, and 4th order 
convergence everywhere for the state vector error 6S. The numerical results are quite close to 
this for a considerable period of time. For example, figures 2 and 3 show typical convergence 
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Convergence of C at t=1000m 
Angular Dependence in +z patch at w=0.12 (r=2.19m) 




FIG. 2 This figure shows the convergence of the energy constraint C for the +z patch of the 
= 0.12 (r = 2.19m) sheh of grid points at t = 1000m. The horizontal axes show the angular 
coordinates in degrees (so (0,0) is the z axis and (45,45) is the "triple point" where the 3 patches 
meet in figure 1). The z axis shows the energy constraint (scaled by the 4th power of the resolution) 
for the 33k-wrmax4, 50k-wrmax4, and 66k-wrmax4 models. 

results for the energy constraint C dXt = 1000m in the 33k-wrmax4, 50k-wrmax4, and 66k- 
wrmax4 evolutions. The best-fitting convergence exponents are about 3.9 for the patch 
interiors and 2.9 for the interpatch boundaries. 

Figures 4 and 5 show typical convergence results for the state vector error 5S at the same 
time in these same evolutions. The convergence isn't as good as for C, but it's still between 
3rd and 4th order everywhere in the grid; the overall best-fitting convergence exponent is 
about 3.8. 



Notice that at larger, C falls off considerably more slowly (~ r~^) than Cabs.cxact (which falls off ~ r"**), 
i.e. the relative constraint violation C/Cabs, exact grows ^ r^. This is due to finite differencing errors 
involving the factors in the (r, p, cr) metric components. It could probably be cured (greatly improving 
the code's accuracy at large r) by factoring those factors out analytically, but I haven't tried this. 
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Convergence of C at t=1000m 
Radial Dependence of r=constant Angular Norms 
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FIG. 3 This figure shows the convergence of various r = constant angular RMS-norms of the energy 
constraint C at t = 1000m. The lower [upper] curves show angular RMS-norms of the energy 
constraint (scaled by the 4th [3rd] power of the resolution) over the patch-interior [interpatch- 
boundary] grid points in each r = constant shell; in all cases oo-norms are within an order of 
magnitude of the RMS-norms shown. The diagonal dashed line shows | [Cabs, exact ||- The two 
vertical lines show the horizon position and the w = 0.12 (r = 2.19m) position used for the angular 
convergence plots. 

F. Time Dependence 

Figure 6 shows the convergence of the energy constraint C and the state vector error SS, 
at the fixed radial position w = 0.12 (r = 2.19m), as a function of time. After an initial 
transient lasting about 250m, the evolution appears to settle down to an almost- stationary 
state (apart from the error waves described below) for some time. Both C and 6S generally 
grow with time, but remain quite small in magnitude (^ 1) throughout the first t = 1000m 
of the evolution. {6S seems to grow at a faster rate than C; this seems to be a gauge 
instability (Alcubierre and Masso (1998); Alcubierre and Schutz (1996).) 

Notice that although the energy constraint C is substantially larger at the interpatch 
boundaries than in the patch interiors (this is visible in figure 2), this difference doesn't 
grow with time, in fact it even decreases a bit. This suggests that the the interpatch 
interpolation is not introducing any instability into the evolution. 

A major exception to the "roughly stationary state" in figure 6 is the sharp spikes in 
C at roughly t = 440m, 680m, 880m, and succeeding times. Figure 7 shows that each 
"spike" is actually the passage of an error wave which originates at the outer boundary and 
propagates inwards at approximately the speed of light. (The error waves are spaced roughly 
250m apart, matching the speed-of-light propagation time from the outer boundary inwards 
to the black hole.) At late times there are also visible error waves propagating outwards 
from the strong-field region, and eventually high-frequency oscillations in C appear in the 
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Convergence of 5S at t=1 000m 
Angular Dependence in +z patch at w=0.12 (r=2.19m) 




FIG. 4 This figure shows the convergence of the state vector error 5S for the +z patch of the 
w = 0.12 (r = 2.19m) shell of grid points at t = 1000m. The horizontal axes show the angular 
coordinates in degrees (so (0, 0) is the z axis and (45, 45) is the "triple point" where the 3 patches 
meet in figure 1). The z axis shows the state vector error (scaled by the 4th power of the resolution) 
for the 33k-wrmax4, 50k-wrmax4, and 66k-wrmax4 models. 

outer part of the grid. 

G. Overlapping versus Just- Touching Patches 

As discussed in section II. D, in my numerical scheme adjacent patches' nominal grids 
must either overlap, just touch, or be separated by a gap of one grid spacing. 

Figure 8 shows a comparison of results from models with adjacent patches just touch- 
ing, versus models with adjacent patches overlapping by ±3 grid points about their mean 
nominal-grid-boundary position (i.e. models with adjacent patches having 7 grid points in 
common in their perpendicular-to-the-boundary direction). The just-touching models are 



The results for the models with adjacent patches overlapped by a fixed angular distance of ±9° are 
generally similar to those for the models with adjacent patches overlapped by ±3 grid points. 
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Convergence of 5S at t=100m and t=1000m 
Radial Dependence of r=constant Angular Norms 
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FIG. 5 This figure shows the convergence of various r = constant angular RMS-norms of the state 
vector error 6S at t = 100m and t = 1000m for the 33k-wrmax4, 50k-wrmax4, and 66k-wrmax4 
models. The curves show angular RMS-norms of the state vector error (scaled by the 4th power 
of the resolution) over all grid points in each r = constant shell; in all cases oo-norms are within 
a factor of 3 of the RMS-norms shown. The diagonal dashed line shows || Cabs, exact ||- The two 
vertical lines show the horizon position and the w = 0.12 (r = 2.19m) position used for the angular 
convergence plots. 

considerably more stable (strictly speaking, less unstable) than the overlapping-patch mod- 
els. In particular, notice that (apart from error waves similar to those discussed in the 
previous section) the just-touching models display roughly constant C for roughly 500m, 
while the overlapping-patch models display steadily rising C in the patch interiors over the 
same time span. 

Another interesting pattern visible in table I and figure 8, is that among the overlapping- 
patch models, the higher-resolution models are more stable (less unstable). It's hard to draw 
any conclusions about a continuum limit, though, because while the just-touching models 
show good 4th (3rd) order convergence in the patch interiors (interpatch boundaries), even 
the highest-resolution overlapping-patch models aren't yet in the asymptotic convergence 
regime. 



The error wave spacing again roughly matches the speed-of-light propagation time from the outer boundary 
inwards to the black hole. 

The formal convergence exponents for the highest-resolution pair of overlapping-patch models are 
about 5(4.5) in the patch interiors (interpatch boundaries). Both of these exponents are greater than 4 
(the order of the finite differencing scheme), so the models can't be in the asymptotic convergence regime. 
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Convergence of C and 5S at w=0.12 (r=2.19m) 
Time Dependence of r=constant Angular Norms 
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FIG. 6 This figure shows the convergence of various r = constant angular RMS-norms of the 
energy constraint C and the state vector error 6S at the fixed radial position w = 0.12 (r = 2.19m, 
as a function of time, for the 33k-wrmax4, 50k-wrmax4, and 66k-wrmax4 models. The points 
show angular RMS-norms of the state vector error (scaled by the 4th power of the resolution) 
over all grid points in each r = constant shell. The thick [thin] lines show angular RMS-norms 
of the energy constraint (scaled by the 4th [3rd] power of the resolution) over all patch-interior 
[interpatch-boundary] grid points in each r = constant shell. 

V. DISCUSSION AND CONCLUSIONS 

Because excision in Cartesian grids is difficult, and polar spherical grids suffer from 
z axis coordinate singularities, there is growing interest in multiple-patch black hole excision 
schemes which allow a smooth r = constant excision surface. As well as greatly simplifying 
excision, such schemes also provide a smooth outer boundary, and (by placing the r = 
constant shells of grid points nonuniformly in radius) they can give high resolution close to 
the black hole while still having the outer boundary relatively far away. 

In this paper I present a detailed description of what I believe is the ffist multiple-patch 
excision scheme for the full nonlinear Einstein equations in 3 -|- 1 dimensions, together with 
sample numerical results from a prototype implementation of this scheme. I use the BSSN 
form of the 3 -|- 1 equations. 

I use a "ghost zone" technique to handle the interpatch boundaries, with all the dy- 
namical fields interpolated into each patch's ghost zones from adjacent patches. I use an 
"inflated-cube" 6-patch system of the form {r x (6 angular patches covering S*^)}. By suit- 
ably choosing the angular coordinates, this allows adjacent patches to always share the 
angular coordinate perpendicular to their mutual boundary, so the interpatch interpolation 
need only be done in the (angular) dimension parallel to the boundary. 

I use each patch's local coordinates to define its tensor basis, so this basis is different 
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FIG. 7 This figure shows 4 frames from a movie showing the time evolution of angular RMS- 
norms of the energy constraint C for the 50k-wrmax4 model. The solid lines (points) show angular 
RMS-norms over the patch-interior (interpatch-boundary) grid points in each r = constant shell; 
in all cases oo-norms are within an order of magnitude of the RMS-norms shown. In each subplot, 
the diagonal dashed line shows ||Cabs,exact|| and the vertical line shows the horizon position. At 
t = 1000m an inwards-propagating error wave is visible at r 100m; by t = 1050m it has 
propagated in to r «i 50m. At t = 1100m the wave is just falling into the black hole, and at 
t = 1150m the wave has almost completely vanished. The full movie is available from http: 
//www. aei.nipg.de/~jthorn/research/mpe/niovies/, and as supplemental information for this 
article on the Classical and Quantum Gravity web site. 

from one patch to another, and hence the interpatch interpolation must include a change of 
tensor basis. This is straightforward for all the BSSN variables except for the F*. The 
aren't tensors or tensor densities, so their change-of-basis transformation law includes terms 
containing spatial derivatives of the BSSN conformal factor 0. Computing these derivatives 
numerically is difficult because this must be done in the ghost zones. I currently solve this 
problem by using a ghost zone twice as wide for the BSSN conformal factor as for the other 
BSSN field variables. This works, but is somewhat cumbersome; it would be interesting to 
explore the alternative scheme I describe (but have not implemented) where the ghost zones 
can be the same width for all the field variables. 

I have implemented my multiple-patch scheme in a prototype numerical code, using 
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Convergence of C at w=0.12 (r=2.19m) 
Time Dependence of r=constant Angular Norms 
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FIG. 8 This figure shows the convergence of various r = constant angular RMS-norms of the energy 
constraint C at the fixed radial position w = 0.12 (r = 2.19m, as a function of time, for the 33k, 
50k, and 66k models, compared to the 33k-overlap3, 50k-overlap3, and 66k-overlap3 models. The 
thick [thin] lines and large [small] points show angular RMS-norms of the the energy constraint 
(scaled by the 4th [3rd] power of the resolution) over all patch-interior [interpatch-boundary] grid 
points in each r = constant shell. 

4th order finite differencing in space and time, and 5th order Lagrange polynomial interpatch 
interpolation. I use ghost zones which are 4 points wide for the BSSN conformal factor 0, 
and 2 points wide for all the other dynamical field variables. My finite differencing scheme 
is a straightforward generalization to 3-D of the 1-D finite differencing scheme I described 
in Thornburg (1999b). The overall accuracy of the scheme is 4th order in the patch interiors, 
and 3rd order close to the interpatch and radial boundaries (though this could be raised to 
4th order everywhere by adjusting the interpolation and radial extrapolation operators). 

In tests of the evolution of octant-symmetry Kerr initial data with various numerical 
parameters, this scheme performs quite well, with evolutions maintaining good convergence 
and preserving the energy constraint near the black hole to better than 1% for ^ 1000m, 
and lasting up to ^ 1500m before crashing. The degradation of accuracy at later times, 
and the eventual crashes, appear to be due to outer-boundary effects, not to any problems 
inherent to the multiple-patch scheme. 

Two obvious areas in which to extend these results would be to drop the octant-symmetry 
assumption, and to evolve more general (non-stationary) initial data. I hope to experiment 
in these directions soon. 

The numerical code I have used for these results is a standalone one, not integrated into 
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1. Radially extrapolate the grid function(s) in the nominal grid to compute their values in the 
nominal-angular-grid part of each patch's inner and outer (radial) ghost zones. 

2. Apply the appropriate symmetry mappings to compute grid function values in the non- 
angular-corner nominal-radial-grid part of each symmetry ghost zone. 

3. Use interpatch interpolation to compute grid function values in the nominal-radial-grid part 
of each interpatch ghost zone. 

4. Do a change-of-basis transformation on each interpatch-interpolated value. (See section II. E 
for details of this.) 

5. Apply the appropriate symmetry mapping to compute grid function values in the angular 
corners of the nominal-radial-grid part of each symmetry ghost zone. 

6. Radially extrapolate the grid function(s) in the nominal grid to compute their values in the 
angular-ghost-zones part of each patch's inner and outer (radial) ghost zones. 

FIG. 9 This figure shows my algorithm for synchronizing ghost zones. 

any standard numerical-relativity toolkit such as Cactus (Goodale et al. (2003)). Such an 
integration will pose some software-engineering problems, but should have large benefits for 
furthering the practical use of multiple-patch schemes. 

How could my multiple-patch scheme be extended to handle multiple (moving) black 
holes? The easiest way appears to be to use an inflated-cube set of patches around (and 
kept dynamically centered on) each (excised) black hole, a Cartesian patch or patches fill- 
ing the space between the black holes, and probably an outer set of inflated-cube patches 
outside both black holes. ^^'^^ This scheme requires multidimensional interpatch interpola- 
tion between the inflated-cube and Cartesian patches, but the experience of Calabrese and 
Neilsen (2004b), as well as the extensive experience with similar schemes in computational 
fluid dynamics (see references cited in section I), suggests that this probably doesn't cause 
any significant problems. Such schemes should be a fruitful area for further research. 

APPENDIX A: Details of the Ghost-Zone Synchronization Algorithm 

In this appendix I give a detailed description of my ghost-zone synchronization algorithm. 

Figure 9 shows the major steps of the algorithm. The main complexity of the algorithm 
(steps 2 to 5) is in the handling of the angular ghost zones. Because this is done indepen- 
dently in each r = constant angular shell of grid points, for the rest of this appendix I drop 
the r coordinate and just describe the angular computations. 

Consider an interpatch boundary between patches p and q, and define (p, a) angular 
coordinates relative to this boundary, ± perpendicular to the boundary and (in each patch) 



Using a Cartesian patch or patches between the black holes avoids the r = singularities of the outer set 
of inflated-cube patches. 

Such a scheme is essentially a 3-D inflated-cube analog of the axisymmetric multiple-patch scheme I have 
previously used for constructing 2-black-hole conformally-flat initial data (Thornburg (1985, 1987)). 
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II parallel to the boundary. (Recall from section TLB that the patches share a common 
_L coordinate, while each has its own || coordinate.) Corresponding to (_L, ||), I define 
(iperp, ipar) as integer grid coordinates to describe the synchronization algorithm, and 
also as integer array indices for subscripting grid functions. 

Since these coordinates are patch-dependent, I use a notation inspired by the C++ 
programming language, where (p. iperp, p. ipar) refers to the (iperp, ipar) coordinates of 
patch p relative to this boundary. Figure 10 illustrates the geometry of the patches, ghost 
zones, and coordinates. 

Suppose p's nominal grid is the rectangle 

(^P®^P' P-ipS-^) £ [P" ^P®^Pnominal,min' P" ^P®^Pnominal,max] ^ [P' ^P^^ nominal, min' P" ^P^^ nominal, max] 

(Al) 

and similarly for q. Consider one of p's (angular) ghost zones gz, which overlaps q's nominal 
grid (that is, patch-p grid function values in gz arc to be obtained by interpolating from 
patch q). Without loss of generality suppose that the coordinates are oriented such that gz 
lies on p's maximum-iperp boundary, so gz's extent in the iperp direction is iperp G 
[p.iperp„„^i„^i^^^ + 1, P-iperp„„^i„^i „,^, + w], where w is the ghost zone width. 

I define gz's extent in the p. ipar direction to depend on the type (symmetry or interpatch) 
of p's adjacent ghost zones on gz's minimum-p.ipar and maximum-p.ipar ends: 

• If p's adjacent ghost zone on an end is a symmetry ghost zone, then I define gz to 

not include their mutual corner. That is, for example, if p's maximum-p.ipar ghost 
zone is a symmetry ghost zone, I define gz's maximum-p.ipar extent to be p. ipar < 
p.iparj^jjj^jj^g^i Part (a) of figure 10 illustrates this case. 

• If p's adjacent ghost zone on an end is an interpatch ghost zone, then I define gz 
to include up to the diagonal line of their mutual corner, arbitrarily breaking ties 
by including the diagonal line itself into the p ghost zone. That is, for example, if 
p's maximum-p.ipar ghost zone is an interpatch ghost zone, I define gz's maximum- 
p.ipar extent at each iperp to be 

^ . / . . X f if ± is p's p coordinate 

p.ipar < p.ipar_,,,_+(iperp-p.iperp_,,,_)-| ^ ^ ^ ^^^^^^^^^ 

(A2) 

Part (b) of figure 10 illustrates this case, with the diagonal line shown dashed. 

Next, I define gz's "patch interpolation region" R in q: This is the set of q grid points 
from which gz's interpatch interpolation will use data. Since the interpolation is only done 
in 1-D along q.ipar lines (these run vertically in figure 10), R has the same extent in the 
iperp direction (horizontally in figure 10) as gz. I define R's extent in the q.ipar direction 
to depend on the type (symmetry or interpatch) of q's adjacent ghost zones on each of q's 
minimum- ipar and maximum-ipar ends: 

• If q's adjacent ghost zone on an end is a symmetry ghost zone, then I allow gz's 
interpolation to use data from that symmetry ghost zone (as well as from q's nominal 



w = 2 for all the numerical results presented here, and for the example in figure 10. 
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grid).^'^ That is, for example, if q's maximmn-q.ipar ghost zone is a symmetry ghost 
zone, I define R's maximum-q.ipar extent to be q.ipar < q.iparj^^jjjj^^^j ^g^^+vj. So long 
as the ghost zone is wide enough, this allows gz's interpolations to be kept centered 
even at gz points within an interpolation-molecule-radius of a symmetry boundary. 
Part (a) of figure 10 illustrates this case. 

• If q's adjacent ghost zone on an end is an interpatch ghost zone (this can only happen at 
a "triple corner" where 3 patches meet), then at present I only allow gz's interpolation 
to use data from q's nominal grid on that end. That is, for example, if q's maximum- 
q.ipar ghost zone is an interpatch ghost zone, 1 define R's maximum-q.ipar extent to 
be q.ipar < q-iparj^^j^^jj^g^j j^^^. Part (b) of figure 10 illustrates this case. Notice that 
for a 4-point or smaller interpatch interpolation molecule, the interpatch interpolation 
can always be kept centered, while for a 5-point or larger molecule, the interpolation 
must be off-centered for points close to the triple corner. 

If interpatch ghost zones are interpolated in some sequential order in step 3 of the 
synchronization algorithm of figure 9, it would be possible to keep more of the inter- 
polations centered by expanding R to include points of the adjacent interpatch ghost 
zone if that ghost zone came earlier than gz in the sequential order, and thus had 
already been interpolated. I haven't investigated this possibility in detail, but (by 
keeping more of the interpolations centered) it might well lead to better numerical 
stability in the time evolution. 



APPENDIX B: Transformation Properties of the 

In this appendix I outline the derivation of the interpatch coordinate-transformation 
law (7) for the BSSN conformal Christoffel symbols f\ 

The transformation law for the non-conformal Christoffel symbols P* is well-known (Mis- 
ner et al. (1973, equation (10.26))): 

P"(p) = X^r(q) + X\Y\,g'^\^) (Bl) 

Substituting the identity (valid so long as we're in a 3-D slice) 

pa ^ g4</>pa ^ 2g''^db(t) (B2) 

into each side of this, then using the (j) transformation law (5) and simplifying, gives the 
final P transformation law (7). 

APPENDIX C: Outer Boundary Conditions 

Basically, I use Sommerfeld outer boundary conditions, applied to each coordinate com- 
ponent of the field variables. However, there's a complication: because the (r, p, cr) tensor 



Since gz is contained in q's nominal grid in the iperp direction, these symmetry-ghost-zone points have 
already been synchronized in step 2 of the synchronization algorithm of figure 9, so it's ok to use their 
values here (in step 3). 
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FIG. 10 This figure shows details of the interpatch interpolation process for ghost-zone synchro- 
nization. The interpatch boundary between patches p and q is shown as a heavy vertical line in 
each part of the figure. In both parts of the figure, iperp increases to the right and ipar in- 
creases downwards, though these coordinates aren't orthogonal. Part (a) shows the case where p's 
adjacent ghost zone on the maximum-ipar end is a symmetry ghost zone (for this example, one 
corresponding to a reflection symmetry across the iperp axis). Sample output points from various 
steps of the ghost-zone synchronization algorithm are labelled. Part (b) shows the case where p's 
adjacent ghost zone on the maximum-ipar end is an interpatch ghost zone. The arrows show 
symmetry operations, the boxed points arc interpatch interpolation inputs, and the circled point 
is an interpatch interpolation output. For the case shown here (a 4-point interpatch interpolation 
molecule), the interpolation can be kept centered, but for a larger molecule it must be off-centered. 
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basis isn't asymptotically Cartesian (for example, asymptotically there are factors in the 
angular metric components), I need to transform to the locally Cartesian basis {x*} be- 
fore applying the usual Sommerfeld condition. In this appendix I use an overbar to denote 
components in a locally Cartesian basis, as in /. 

For each field variable /, I actually apply a Sommerfeld condition not to / itself, but 
rather to the deviation 6f = f — /background of / from its (time-independent) value in a 
specified background slice. The Sommerfield condition for / is thus dt{r^f) = —cdr{r"'6f), 
where c is the propagation speed and n the falloff exponent, or equivalently 

— ' — 11 — 

dtf = -c dr5f + -5f (CI) 
L r J 

For 3-scalars such as a and K I can apply this boundary condition directly. The other 
BSSN variables need to be transformed to locally Cartesian coordinates for the Sommerfeld 
condition (CI) to be apphed. Proceeding analogously to section II. E, I define the transfor- 
mation matrices 

X\ = ^ (C2a) 

y\ - ^ (C2b) 

Y\b = -^^^ (C2c) 

These are straightforward to compute analytically from the coordinate definitions (1). Note 
that the X°'i and Y'^a matrices are inverses. 
The BSSN conformal factor transforms as 

4> ^ 4^ - \\og\Y\ (C3) 

where Y = det [?'a] • Because Y = Y(x"') is a known coefficient, not a dynamical variable, 
S(l> = 6(f), so for purposes of the Sommerfeld condition (CI) I can treat as a 3-scalar, and 
apply (CI) directly. 

The BSSN gij and Aij transform as 

4- = \Y\'/'X\X',Sab (C4) 

where Sij G {gij.Aij}. Substituting this into both sides of the Sommerfeld condition (CI), 
multiplying both sides by \Y\~'^^^Y\Y^ d and summing over i and j, and simplifying, then 
gives the boundary condition 

dtSij = -c 



drSSij + -SSij + Z\5Skj + Z^^SSk, + l{dr log \Y\)5Sij (C5) 



where Z% = Y%drX^i. 

The BSSN transform in a more comphcated way: Applying the interpatch transforma- 
tion law (Bl) to the coordinates x"(p) and a;*(q), and using the identity (B2) and its inverse 
in each of these coordinate systems, gives the transformation law 



|y|-2/3 



yij^a _ 2Y\g-%<f> + 2Y\Y\r%<P - Y'aW'^ + \Y\Y\{d, log \Y\)r 



(C6) 
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Substituting this into both sides of the Sommerfeld condition (CI), multiplying both sides 
by and summing over i, and simplifying, then gives the boundary condition 



+ X%Y\^dtr' + \y\{dj log \Y\)dt~g^ 



be 



-2 {drSg'nfl) + Sg'^drbcl) + drg'^ScI) + g'^drbScj)) 
+ 2W\Y\ {Sr%<P + r%S(f>) + 2{drY\) {5~g"'dj(f) + g''dj5(f)) 
+ 2Y\ {dr5~g^''dj(t) + 5~g^''drj(t> + drg^''dj5(t) + g'"'drj5(f)) 



- lW\Y\idj log \Y\)5r' - ldrY\{dj log \Y\)5~g' 

- lY^bidrj log 1^1)5^"^ - lY\{dr log \Y\)dr5g' 



~bc 



be 



+ 



(- 



ldAog\Y\ 



6r^ - 2 {6g^^db(/) + g^''dbS(j)) + 2Y\ {Sg^^dj(j) + g^^d^ 



- X\Y\bbcf' - ^^\{di log \Y\)8~g 



be 



(C7) 



where W"'b = X"'idrY^b and W°'bc = X"'id,Y'^bc- The right-hand-side terms involving time 
derivatives can be evaluated by using the other boundary conditions: Differentiating the 
Sommerfeld condition (CI) for gives 



r 



+ 



—n/r^ 




if i = r 
otherwise 



(C8) 



while dtg"'^ can be evaluated by applying the matrix identity dg"'^ — —g°''^g^^dgij to the 
boundary condition (C5). 

For a, and 0, I use a propagation speed c given by the Bona-Masso lapse propagation 
speed (9), while for g^j, A^j, and I use a propagation speed given by the outgoing light- 
cone speed (xj-^fg^ — (3'^ . I use a falloff power of n = 1 for all the Sommerfeld boundary 
conditions. 

For historical reasons, all the numerical results reported here were obtained using a simple 
Dirichlet outer boundary condition SjF* = instead of the Sommerfeld condition (C7). This 
probably contributes to the outer boundary instabilities I see, and I plan to switch to the 
Sommerfeld condition for future computations. 
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